The large system asymptotics of persistent currents in mesoscopic quantum rings 
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We consider a one-dimensional mesoscopic quantum ring filled with spinless electrons and threaded 
by a magnetic flux, which carries a persistent current at zero temperature. The interplay of Coulomb 
interactions and a single on-site impurity yields a non-trivial dependence of the persistent current on 
the size of the ring. We determine numerically the asymptotic power law for systems up to 32 000 
sites for various impurity strengths and compare with predictions from Bethe Ansatz solutions 
combined with Bosonization. The numerical results are obtained using an improved functional 
renormalization group (fRG) method. We apply the density matrix renormalization group (DMRG) 
and exact diagonalization methods to benchmark the fRG calculations. We use DMRG to study 
the persistent current at low electron concentrations in order to extend the validity of our results 
to quasi-continuous systems. We briefly comment on the quality of calculated fRG ground state 
energies by comparison with exact DMRG data. 
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I. INTRODUCTION 

Quantization of the magnetic flux was first observed 
in superconducting cylinders [l], [2] ■ In normal electronic 
systems the flux quantum 0o = hc/e determines the pe- 
riod of the persistent currents, which can be observed in 
mesoscopic one-dimensional disordered metallic rings Q . 
Persistent currents in metallic rings are induced only if 
the circumference of the ring is not larger than the elec- 
tron phase coherence length and/or the electron's mean- 
free path. Such conditions can be easily obtained for 
sufhciently low temperatures at which inelastic scatter- 
ing from phonons is suppressed. Normally, this happens 
at temperatures below 1 K and ring sizes of about 1 /im. 
Impurities (disorder) in such metallic rings act as elastic 
scatterers, and it is not surprising that persistent currents 
can flow despite the non-zero resistance of the rings. 
Interest in persistent current Tjhenomena remains 
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strong both experimentally 
theoretically [ll [H, [H, [Is 

[2^ . [25I . The main reason for this is the fact that there 
still is a discrepancy between experiment and theory. Ex- 
periments [1, 11] show persistent currents which are two 
orders of magnitude larger than theoretical predictions 
for rings filled with non- interacting electrons [H, [l^, . 
If electron-electron interactions are included, persistent 
currents increase for repulsive as well as attractive elec- 
tron interactions 0,113 113' however, the calculated cur- 
rents were still about five times smaller than the experi- 
mental data. Recently, Bary-Soroker et al. [11] provided 
an explanation for the magnitude of the persistent cur- 
rents if magnetic impurities and the attractive interac- 
tions are considered. 

Considerable theoretical effort has been invested into 
the question whether the interplay between strong 
electron-electron interactions, electron density, and dis- 
order strength can explain the large persistent currents 




observed experimentally. Disorder usually gives rise to 
elastic scattering and consequently should reduce the 
persistent current. Bouzerar et al. (l5| claim that both 
electron-electron interactions and disorder decrease the 
persistent current, whereas the self-consistent Hartree- 
Fock calculations by Cohen et al. [11] suggest that the 
current may be enhanced by a moderate disorder due to 
screening effects. Other theoretical studies conclude that 
electron-electron interaction could reduce or enhance the 
current de p endin g on the type of disorder present in the 
system [lj,[ll,[l3]- Abraham and Berkovits [13] conjec- 
tured that the increase would be rather small and one- 
dimensional spinless models at half-filling cannot explain 
the results of single ring experiments. 

Among the theoretical methods used to study per- 
sistent current phenomena arc Hartree-Fock approxima- 
tions [r^ [Tsj . configuration interaction and quantum 
Monte Carlo simulation s Isol . Isij ] , Bosonization [3] , con- 
formal field theory [l9l 1201 ] , the functional renormaliza- 
tion group (fRG) |21| . and the density matrix renormal- 
ization group (DMRG) [2ll,[2i,[2|. 

In the present paper we study the interplay between 
electron-electron interactions and a single impurity in a 
one-dimensional system. For one-dimensional interact- 
ing systems, the presence of a single irnmirity is known 
to affect their physical properties [H, jH, [SJ, [s^. We 
restrict our investigations to the Luttinger liquid phase 
which is characterized by a power law decay of correlation 
functions in a sufficiently large system [36] . Our aim is to 
investigate the asymptotic behavior of the persistent cur- 
rent / as a function of the ring length L, from which the 
power law exponent may be obtained / oc L^"^^^. The 
exponent as = — 1 is related to the Luttinger liquid 
parameter K [H, [s^, [s^ and depends on the electron- 
electron interaction, but it does not depend on the im- 
purity strength if L — > cxd. 

It is known that the asymptotics of the currents is 
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reached for smaller system sizes if sufficiently strong im- 
purities are considered [H, [13, HI]. Nevertheless, it is 
necessary to investigate relatively large systems in or- 
der to reach the asymptotics. We therefore need suit- 
able many-body techniques which enable such calcula- 
tions. We have chosen the functional renormalization 
group [UliO, EH 113, m, which proved to be a rather 
useful tool in a previous study [21[. However, for nu- 
merical reasons the asymptotic regime could not really 
be reached in that work. In the present paper we im- 
prove the fRG technology such that system sizes up to 
about 32000 lattice sites can be investigated. In order 
to benchmark the fRG results we use DMRG calcula- 
tions [ii,!!!,!!!,!!^. 

We also study the asymptotic power law decay of the 
persistent current at low electron concentrations in order 
to clarify differences between discrete lattice models (in- 
vestigated in this paper) and continuous systems studied 
by quantum Monte Carlo and configuration interaction 
methods [s^, [3l|. For this reason, we apply DMRG to 
our model away from half-filling (n = i), in particular, 
we study the case of quarter-filling (n = i) as well as 
n — ^ and n — j^. 

The paper is organized as follows. In Sec. 2 we de- 
fine the model for our calculations and briefiy review the 
theory of the persistent current. We introduce the fRG 
method in Sec. 3 and related Appendices. Section 4 is 
devoted to the DMRG method with a short discussion of 
its limitations. Results are presented in Sec. 5 gathering 
the numerical data obtained by fRG, DMRG, and exact 
diagonalization. We calculate the persistent current and 
extract the effective exponent acs related to the power 
law decrease with respect to the ring size L. Finally we 
show DMRG results for systems with low-filling. The 
results are summarized in Sec. 6. 



II. MODEL FOR SPINLESS FERMIONS 

We consider a quantum ring of interacting spinless elec- 
trons at zero temperature. The ring is pierced by a mag- 
netic fiux (f> and contains one single impurity. This system 
is modelled by the tight-binding Hamiltonian 

L L 

n = -tY^ (clc^+ie""^/-^ + h.c.) +U^ntni,+i + V ni , 

£=1 f.=l 

written in terms of the electron creation and annihila- 
tion operators c\ and q as well as the electron density 
operator ni — c\ci. 

The Peierls factor e~^'^^^ with = 27r$/$o describes 
the influence of the magnetic field in terms of the mag- 
netic flux $ [4^. The flux quantum $0 = ^c/e is set to 
one in the following. The hopping amplitude t between 
neighboring sites is also set to one in order to define the 
energy scale. The nearest-neighbor Coulomb interaction 
is denoted by [/. An (external) on-site impurity potential 
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FIG. 1: (Color online) The five lowest-laying single-particle 
energy levels for a non-interacting system as functions of the 
magnetic flux for a ring of size L — 128. The three panels 
correspond to V = 0, V = 0.1, and V = 1.0, respectively. 



with strength V is placed on the lattice site £ — I along 
the ring {£ — 1,2,...,L). Periodic boundary conditions 
are imposed by i + 1 = 1. 

Without electron-electron interaction {U = 0) and im- 
purity {V = 0) the single particle energy levels cal- 
culated from the Hamiltonian Eq. ([T]) depend quadrat- 
ically on the magnetic fiux (f). The persistent current 
Im = —Vm/L at energy level m with the wave vector 
k — 2tt {(f) + m) /L is calculated from the electron veloc- 
ity V,n, 



Vmi4>) = 2tt 



dk{(f>) 



L- 



(2) 



Summation over all occupied energy levels then yields 
the total persistent current at temperature T = 0, 



E 



i9e„(0) 



(3) 



The upper panel in Fig. [T] shows the parabolic energy 
band structure {V — 0, U ^ 0). Consequently, the per- 
sistent current is proportional to the magnetic fiux </) and 
inversely proportional to the length L of the ring. For an 
odd number of electrons in the ring one finds 



VF 



-TT < < TT , 



(4) 
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whereas for an even filling the current is given by 



+ 71"; 



- TT; 



-TT < (/)< 0, 

< (h<n 



(5) 



with v-p being the Fermi velocity. 

An impurity in the ring rounds off the cusps seen in 
this function as shown in the central and lower panels of 
Fig. [1] for a weak and an intermediate impurity, respec- 
tively. For strong impurities and odd filling, one finds 
/ oc — sin(0)/7rL [ig. The couphng strength V of an 
on-site impurity can be related to the physically more rel- 
evant transmission coefficient at the Fermi wave vector 



r(fcF) = 



4^2 ^ Y2 



(6) 



which holds at half-filling in the non-interacting 
limit [1^. With this relation, the strengths of various 
impurities V used in this paper may be compared. 

If the electron-electron interaction U is switched on, 
the many-body energy levels are shown in Fig. [2] for 
U = 1 and various impurity strengths V. The persis- 
tent currents are obtained from the relation 
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FIG. 2: (Color online) The five lowest-laying many-body en- 
ergy levels obtained by DMRG as a function of the magnetic 
flux (p for a ring of size L = 128 at half-filling and electron- 
electron interaction U = 1. Upper panel: V = 0, center panel: 
V = 1, lower panel: V — 2. 
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where Eo{4>) is the many-body ground-state energy at 
zero temperature as plotted in Fig. [21 In general the 
system is a Luttinger liquid. However, at half-filling it is 
in the Luttinger fiquid phase only for |[/| < 2. At [/ = 2, 
the system exhibits a charge-density wave instability, and 
phase separation characterizes the system for U < —2. 

In this paper we focus on the Luttinger liquid phase 
for which Bosonization techniques and additional approx- 
imations predict that asymptotically (i.e. for large L) 
the persistent current decays algebraically with increas- 
ing system size. 



I ocL- 



(8) 



The exponent ae is a function of the electron-electron in- 
teraction (regardless of the impurity strength V). Only 
at half-filling, the exponent ae can be obtained analyti- 
cally from a Bethe Ansatz solution [U [3^, [s^l 



ttB 



2 / U 
— arccos 

TT V 2 



(9) 



Throughout this paper we often consider the case U — 
1 with the corresponding ae — ^. We determine as 
from the Hamiltonian Eq. ([1]) using three different many- 
body techniques: fRG, DMRG, and exact diagonalization 
(ED). 



III. FUNCTIONAL RG SCHEME FOR 
INTERACTING FERMIONS 

A scheme for one-dimensional Fermionic systems based 
on the functional renormalization group, has been devel- 
oped by Meden, Metzner, Schonhammer, and collabora- 
tors fsil based on seminal work by Wetterich [i^ and 
Morris [4l| . It has been applied to study transport prop- 
erties of quantum wires [5l|, [52, [SJ] and rings |21j] . 

The fRG scheme of the present work is similar to the 
one applied in the references cited above, but uses a 
somewhat different truncation procedure. Details of our 
method have been presented recently [s^ . Here we re- 
view the essential steps relevant for a system of spinless 
interacting Fermions. 

The effective average action Ffc for interacting 
Fermions evolves according to the fRG equation [40, ^ 



-iTr 
2 



Rk 



dk 



(10) 



The system of L spinless electrons is described by an 
L-component vector of Grassmann variables (/>(t) = 
((/)i(t), . . . , 0i(T)), which describes the evolution of the 
interacting Fermionic system in imaginary time r. The 
functional derivative of F^ with respect to 4>* and is de- 

(2) 

noted by . The trace in Eq. (|10p is to be performed 
over the Fermionic states j = 1, . . . , L. The regulator 
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Rk is introduced in order to suppress thermal and quan- 
tum fluctuations at energy or momentum scales k larger 
than any physical scale relevant for our problem. With 
decreasing k, the regulator gradually "switches on" such 
fluctuations until they are fully included at A: = 0, i.e. 
at fc = the regulator Rk vanishes. As a regulator we 



choose Rk 



■ w 



') with C a large constant and 



it satisfies all requirements for a useful regulator as dis- 
cussed in more detail in Ref. js^l- It has an additional 
advantage that the integration over the frequencies w can 
be done analytically. 

In order to solve the functional differential equa- 
tion (jlOp a particular functional form for Ffc in terms 
of the Grassmann variables , is assumed 

Ffc W dr ^ C (r) (r) ^U^ (0* (r) , 0(r)) , 

(11) 

where the 'effective potential' lA^ does not depend on 
derivatives of the Grassmann variables with respect to 
the imaginary time r, i.e. it is represented as a Grass- 
mann polynomial in terms of 0* and <i) 



L 



Uki4>*,4>) = ak,o + ^ -I- afcjj+i(/)*(/)j 



-1-1 



+ 0.k.j,j-l4'j(f'j-l 



L 



(12) 



with the understanding that 00 = (pL and (/>l+i = 01- 
Of course, in general this polynomial should have higher 
order terms, but those are neglected in the hope they are 
small. 

For half-filling, the density-density interaction strength 
Uk renormalizes according to 



Uk = 



u 



1 



2tt 



(13) 



as was suggested in Ref. [SOj. This result can be eas- 
ily derived within the fRG scheme applied here by using 
analogous approximations for the two-particle vertex as 
in Ref. [soj . (The interaction parameter U is given in the 
Hamiltonian Eq. ([1])). 

Inserting this Ansatz for the effective average action on 
both sides of the flow equation, performing the integra- 
tion over the frequencies lu, and comparing terms with 
the same Grassmann structure, it is straightforward to 
obtain the flow equations for the coefficients ak,o and 



k,jj 



L^,.«M„,,(lg-(.A)) 

\=±k ^ ^ 



Uk 



si: i: 



•9fc.0+s)0+s)(«A), 



a 



k,]U±l) 



\=±k 



(14) 



with gfc(iA) = (afc + iXl) and — {ak.jj')- The con- 
vergence factor e*'^"^ is only needed at the beginning of 
the flow from infinity down to a large constant fco- 

At fc = oo the effective average action Tk is given by 
the classical action S (iol . [ssj , which is determined by the 
Hamiltonian Eq. Therefore, the initial conditions at 
k = oo for the solution of the flow equations ([T^ are 
directly obtained from Eq. ([T]), 



<^oo j'(j±l) 



= -VSji + fiS 



(15) 



with the understanding that aL,L+i = Q-l,! and oi^o = 
ai,L- For convenience we have added a chemical potential 
term —jini to the Hamiltonian. Our initial conditions 
differ from those of Ref. |2l|: in the present work the 
dependence on the magnetic flux <f> resides in the initial 
conditions for the off-diagonal ajj' only, while in Ref. [2lj 
this dependence partly resides in the free propagator, cor- 
responding to a different definition of the kinetic energy. 
Moreover, here we also include the chemical potential 
into the ajji and not in the free propagator. 

In practice we must start the renormalization flow at 
a finite k = kg. The initial fiow from k = oo to k ^ kg 
is obtained from an analytical solution of the Eqs. (|14p . 
For convenience this calculation is briefly reviewed in Ap- 
pendix]^ with the result 



ako,o 

^kojj 



_ V L fU 
2" 2" V 2" 



(16) 



\=±k s=±l 



At half-filling the chemical potential is given by /i = [/ 
because of particle-hole symmetry. From the ground 
state energy Eo{(j)) = ao{k = 0) -f fi{n), the persistent 
current is calculated using Eq. ([3]). 

For given L, Eq. represents a (large) set of 3L — 1 
non-linear coupled complex-valued ordinary differential 
equations. The calculation of the right-hand side of these 
equations requires the inversion of a (potentially) large 
cyclic tridiagonal matrix as well as the computation of 
the logarithm of the determinant of such a matrix at 
each step of the numerical integration of the differential 
equations. To accomplish this in an efficient manner is 
described in some detail in Appendix IB] 



IV. DMRG 

The density matrix renormalization group is a numer- 
ical technique for the diagonalization of the very large 
matrices typically encountered in quantum many-body 
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calculations. The technique is described in detail in 
Refs [H,!!!, mill]. We use DMRG for the calculation of 
the ground-state energy -Eq (</')• In our case, the matrices 
to be diagonalized are complex- valued due to the Peierls 
factor e^*"^/^ entering Eq. ([T]). Treating such complex- 
valued systems by DMRG does not lead to numerical 
instabilities, apart from a few rare cases which occur at 
very low electron fillings. Here, however, the numerical 
diagonalization routines for complex matrices experience 
convergence difficulties for very large system sizes. Also, 
the superblock diagonalization routines seldom require 
more than 10"* cycles for convergence as compared to the 
standard ~ 10^ cycles typically necessary for real-valued 
matrices. The memory requirements increase by a factor 
of 1.8 and the calculation time raises by a factor of 2.1 
compared to a standard real-valued DMRG. 

We kept the DMRG truncation error as small as 
e < 10~^ for system sizes L > 128, for smaller systems 
e ~ 10~^^ could be achieved. The number of states kept 
are left to vary such that the above truncation error con- 
dition could be satisfied. It is well known that the effi- 
ciency and accuracy of DMRG decreases substantially if 
periodic boundary conditions are imposed. We checked 
the accuracy of our results comparing with data obtained 
from the exact diagonalization of the Hamiltonian for sys- 
tems with up to L = 24 sites at various fillings. 

We would like to point out that the calculations for 
strong impurities > 10 and for large rings with L > 100 
require an extremely careful treatment because differ- 
ences between the ground state energies for different mag- 
netic flux values are extremely small \Eo{(f) = 0) — i?o(0 = 
TT)\/\Eo{(j) = 0)1 < 10"^. For this reason, both the 
number of states kept and the number of the finite-size 
method sweeps must be sufficiently large. 



V. RESULTS 

In this section we present and analyze numerical data 
obtained by three different methods: fRG, DMRG, and 
ED. We study the dependence of the persistent currents 
on the ring size L, the impurity strength V, the magnetic 
flux 0, and on the electron concentration n. Varying the 
electron-electron interaction U within the Luttinger liq- 
uid regime — 2 < [/ < 2 does not change physical prop- 
erties qualitatively as we shall see below. However, the 
weaker the electron-electron interaction the larger is the 
system size L needed to reach the asymptotic power law 
decay of the persistent current [SQi] . 

The ground state energies Eo{(j)) are calculated for < 
(j> < IT and extended to the first Brillouin zone — tt < (f> < 
TT using the reflection symmetry of the energy around the 
origin = (cf. Figs. [T] and [2|). The persistent current in 
Eq. ^ obtained by DMRG and ED is calculated using 
numerical differentiation. The Fourier coefficients Ik of 



the currents 

oo 

/(0) = ^4sin(fc0) (17) 
fc=i 

can be calculated without numerical differentiation using 
an integration by parts 

/fe = — / d(j)Eo{c^)cos{k^). (18) 

Nevertheless, small numerical errors in the calculation of 
the ground state energy may significantly affect the cal- 
culation of the power law exponent ag (cf. Eq. ([5])). The 
Fourier analysis for strong impurity strengths in Eq. p7D 
is well justified due to the sine-like shape of the current, 
and the higher coefficients Ik, k > 2 decay to zero with 
increasing L faster than Ii. However, if the impurity 
y <C 1, the first Fourier coefficient Ii may not suffice 
to characterize the decay of the persistent current accu- 
rately. 

For system sizes L < 24 and for impurities V < 100, 
ED and DMRG yield essentially exact ground state en- 
ergies in mutual agreement. The ground state energy 
calculated from fRG differs from the exact results of ED 
and DMRG at large interactions, as shown in Fig. [3] for 
interactions within the Luttinger liquid regime. This de- 
viation signals the drastic approximations involved in the 
fRG method used here, in particular, the discarding of 
higher order vertex functions. However, the shape of the 
ground state energy as a function of the magnetic flux 
can be reproduced quite well by this method as can be 
seen from Fig. [H Therefore, we can use fRG calculations 
in order to obtain the persistent currents for systems as 
large as L = 3 • 10^. For very strong impurities {V :s> 10) 
DMRG and fRG may run into numerical difficulties and 
we analyzed such cases by ED as is discussed in more 
detail below. 



A. fRG and DMRG analysis at half-filling 

Figure m shows persistent currents as functions of the 
flux calculated at half-filling (n = i) using DMRG and 
fRG for various ring lengths L and for an intermediate 
impurity strength V — 2 (corresponding to transmission 
T = i). The data indicate that for small and interme- 
diate ring lengths L, the DMRG and fRG calculations 
agree rather well. DMRG calculations are hardly feasi- 
ble for L > 256, and only fRG data are available in order 
to study the large L limit. 

Bosonization theory including additional approxima- 
tions predict a power law decay of the persistent cur- 
rent at large L, cf. Eq. ([8]). In Fig. [5] we compare this 
asymptotic decay (dashed lines) with the numerically de- 
termined first Fourier coefficient of the persistent current 
LIi in a logarithmic plot. From Bosonization one expects 
that LIi decays asymptotically as L~"^. The open cir- 
cles represent the fRG data and the filled triangles the 
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FIG. 3: (Color online) Ground state energies at </!> = calcu- 
lated by DMRG (filled triangles) , ED (open circles) , and fRG 
(crosses) as functions of U within the Luttinger liquid regime 
for L = 16 and V = 1. 




FIG. 4: (Color online) Persistent currents calculated by fRG 
(fines) and DMRG (symbols x) at half- filling for f/ = 1, 
V = 2, and various ring lengths L. The currents with the 
largest amplitude correspond to L = 4 and then the ampli- 
tudes decrease with increasing L. 



DMRG results. Again we see the good agreement be- 
tween both methods, and it is obvious that the asymp- 
totic decay does not depend on the impurity strength (if 
at least V > 0.1), which ranges from rather weak V — 0.1 
(T = 0.9975) to very strong V = 100 (T = 0.0004). The 
asymptotic regime is reached at smaller L for stronger 
impurities. 

Notice a tiny deviation of the fRG data from the 
dashed lines with the common power-law exponent q;b 
at large L. This is due to the approximations involved in 
the fRG method and is known from studies of quantum 
wires {2^. However, there are also numerical limitations: 
for large systems the differences between currents calcu- 
lated for different magnetic fluxes are so tiny that they 
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FIG. 5: (Color online) The first Fourier coefficient h versus 
the length of the ring L for U = 1 and various V. The open 
circles and the filled triangles, respectively, correspond to the 
fRG and DMRG data. The dashed straight lines show the 
asymptotics qb = | . The upper a:-axis shows the exact num- 
ber of the sites L at which the calculations were performed. 
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FIG. 6: (Color online) The effective exponent aeff versus the 
inverse of logjQ (L) for the interaction U = 1 and various im- 
purities V. The open and filled symbols correspond to fRG 
and DMRG data, respectively. 



can not be resolved by the floating point data represen- 
tation of the computer, and the current cannot be calcu- 
lated reliably. This numerical limitation also prevents us 
to calculate even larger systems. 

In order to quantify the deviation of the numerical data 
from the expected power law, we define an effective ex- 
ponent 



91ogio(L/i) 



(19) 



91ogio(L) 

which is shown in Fig. [H] as a function of l/log]^Q(i) 
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FIG. 7: (Color online) The persistent current versus L for a 
variety of V's obtained by ED at [/ = 1 and half-filling. 
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FIG. 8: (Color online) Dependence of the effective exponent 
on 1/ logjQ (L) for the same parameters as in Fig. [T] 



The data labeled by the open and filled symbols are ob- 
tained from the fRG and DMRG results shown in Fig. [5] 
using a numerical derivative. If the impurity is weak 
{V = 0.1), the effective exponent aeB deviates signifi- 
cantly from ae = | (the horizontal dashed line) . It is ex- 
pected that substantially larger systems would need to be 
considered to approach the asymptotics in this case. For 
V > 0.5, the numerically determined exponents tend to 
saturate at Qcff ~ 0.35 which is about 5% larger than the 
expected ae = -j- Our result is in good agreement with 
the leading [/-behavior, where the correction to the expo- 
nent for the quantum wire is known to be - = 0.318 SI. 
As already stated, we attribute this deviation to the ap- 
proximations involved in the fRG. Similar results were 
obtained in a study on open chains [53l | , where the same 
power law decay is extracted from the decay of the spec- 
tral weight near the impurity site. 

From our data analysis it emerges that the effective 
exponents for larger L are also non-negligibly infiuenced 
by the numerical limitations discussed above, which are 
seen in the figure as small irregularities in the approach 
of the numerical data to saturation. Notice that we are 
dealing with extremely tiny effects observed in the ex- 
ponent Qfeff- Thus, extremely accurate calculations of 
the ground state energy is the key to the analysis of the 
effective exponent. There is another feature of the effec- 
tive exponent worth mentioning: it shows a characteristic 
minimum at lengths 8 ^ i < 100 within the range of V. 

The effective exponent determined from the DMRG 
data (full symbols) is systematically below the fRG re- 
sults for all but the very small systems. Since the DMRG 
data are expected to be accurate, one may conjecture 
that ttcff approaches ub = for L ^ oo, ii DMRG cal- 
culations for such large systems would be feasible. 



V 

1 2 4 10 10" 10^ lo" lo' lO*" 




FIG. 9: (Color online) Relative difference between the effec- 
tive and asymptotic exponents as a function of the impurity 
strength V for several sizes L. 



B. Non- monotonous behavior of the effective 
exponent 

Another interesting feature can be seen in Fig. [6l 
Both fRG and DMRG data show that the result for the 
strongest impurity V = 100 appears to be out of sequence 
if compared to the other F's. To elucidate this further, 
we performed a series of calculations using ED at smaller 
sizes L. Figure [7] shows the first Fourier coefficient of the 
persistent current Ii as a function of L on a logarithmic 
scale for on-site impurities in the range 1 < < 10^. In 
order to gather the data, the Fourier coefficient is scaled 
with VL in Fig. [T] Furthermore, we plot the absolute 
value of the Fourier coefficient in order to remove the 
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20 24 



FIG. 10: (Color online) Persistent currents as a function of 
the system size L for magnetic flux 4> = ^, various electron 
concentrations n, and interactions U . Plotted is the absolute 
value of the persistent current in order to remove the sign 
change between even and odd electron fillings. 



even-odd effect discussed in Section II. The straight dot- 
dashed hne corresponds to a-Q = \- 

Figure [8] shows corresponding calculations of the effec- 
tive exponent for the same set of impurity strengths. 
The dotted lines only serve as guides to the eye and 
should not be taken as extrapolations. Again, it is ob- 
vious that the data for the cases V > 10^ appear to be 
out of sequence as in Fig. [6l In Fig. [9] we plot the rela- 
tive difference Aq, — (ae — ctcs)/c(B hi order to quantify 
the dependence of the effective exponent on the impurity 
strength V for several system sizes L. A minimum of 
Aq appears at F ~ 4 for J7 = 1. We, therefore, identify 
three regions in the figure (for reachable system sizes): 
in the first region, < y < 4, we observe the stan- 
dard behavior, i.e., the stronger the impurity, the faster 
the convergence of aeff to ae. In the the second region 
which starts at around V ^ A and ends around V — 10"^, 
this behavior is reversed. In the third region V > 10'^ 
the effective exponent acs saturates and does not change 
any more with increasing impurity strength. 



C. DMRG study at low fillings 



10 



10" 



10" 
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FIG. 11: (Color online) Persistent currents versus system size 
L for 17 = 1 at various electron fillings n. The open and 
full symbols correspond to the impurity V — 2 and V = 10, 
respectively. The slopes of the dashed lines coincide with the 
asymptotic exponents as- 
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FIG. 12: (Color online) The effective exponents versus the 
inverse of logj^g(L) for the data displayed in Fig. 1111 The hor- 
izontal dashed lines correspond to the asymptotic exponents 
Qfeff (w) with n as listed in the legend. 



Studies of the persistent current at low fillings try to 
address the question of differences between discrete lat- 
tice models and continuous models of one-dimensional 
systems. Such questions have been addressed in a re- 
cent comprehensive numerical study [13, [3l|. Discrete 
models mimic continuous electron gas models if the elec- 
tron concentration n is such that the cosine dispersion of 
the discrete system can reasonably well approximate the 
parabolic dispersion of the electron gas. Moreover, lat- 



tice models violate Galilean invariance [5^ . In a Galilean 
invariant model without impurity {V = 0) the persistent 
current should not depend on the interaction U. 

First we check whether the interaction dependence of 
the current really disappears at lower fillings: Figure [TOl 
shows persistent currents as a function of system size 
L for U = 0, U — 1, and U — 2 a,t various fillings 
n — ^, J, and i. The currents at half-filling (full sym- 
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bols) strongly depend on the electron-electron interac- 
tion, while the electron gas is better approximated at 
quarter-filling (shaded symbols). A negligibly small de- 
pendence on U is found at n = i (open symbols) where 
Galilean invariance is almost satisfied. 

Figure [Tl] shows the decay of the persistent current 
as a function of the ring size L for electron concentra- 
tions JT- = ^, 3:, |, and on a logarithmic scale for 
V = 2 (open symbols) and y = 10 (full symbols). The 
asymptotic values ae {n) expected from the Bosonization 
analysis [36, 50, 57] arc indicated by the dashed straight 
lines corresponding to exponents aB{n ~ \) — -j, 
aB(i) = 0.18383, aB{\) = 0.08825, a^ijo) = 0.042940. 
On the scale of the figure the data seemingly reach the 
asymptotics. A more detailed picture of how the effective 
power-law exponents aoff approach the expected asymp- 
totic values ttB is shown in Fig. 1121 The upper and lower 
panels show the effective exponent calculated for y = 2 
and V = 10, respectively, at various fillings n. It is ob- 
vious that the acff(ri) do not converge to a^in), and 
substantially larger system sizes L would need to con- 
sidered in order to see convergence. From our data we, 
therefore, cannot support the suggestion [s^, HH that the 
asymptotic regime can be reached at a small number of 
electrons nL in continuous models. 



region V > 10^ the dependence of the effective exponent 
ttoff on the impurity strength V saturates. 

Since discrete lattice models are not Galilean invariant, 
we decreased the electron filling down to n = such 
that the discrete Hamiltonian mimics a continuous elec- 
tron gas with a quadratic dispersion law. From our data 
we see that even at low fillings, i.e. for quasi-continuous 
models, we need large systems in order to reach the ex- 
pected power law exponents. 
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APPENDIX A: SOLUTION OF THE FRG 
EQUATIONS FOR LARGE k 



VI. SUMMARY 

We studied the asymptotic power law decay of per- 
sistent currents using three different numerical methods 
(fRG, DMRG, ED). To this end we improved the fRG 
method so that systems with up to 3 • lO'* sites could be 
treated. We used DMRG and ED in order to benchmark 
the fRG results and confirmed a sufficiently good agree- 
ment between fRG, DMRG, and ED calculations. We 
confirmed that the fRG is a suitable method for extract- 
ing the asymptotic Luttinger power law exponent ae- 

We found that ae which describes the decay of the 
persistent current is about 5% overestimated by the fRG 
results compared to expectations from Bosonization and 
additional approximation. This deviation we attribute 
to the approximation of the fRG method, and the effec- 
tive exponent is known to be correct to the leading order 
in the electron-electron interaction from the calculations 
of the same exponent for quantum wires. The accurate 
DMRG analysis cannot, unfortunately, be extended to 
sufficiently large systems. However, there are indications 
that the expected asymptotic values could be reached, if 
larger systems could be calculated. It is found that the 
effective exponent for a fixed impurity strength V shows 
a typical minimum for lengths 8 ^ L ^ 100 at half-filling. 

Generally, it is known that the stronger the impurity, 
the faster the asymptotic regime is reached. However, 
we identified three regions for U =1: in the first region, 
< y < 4, we observed the standard behavior. (The 
case V = A corresponds to T = i.) In the the second 
region 4 < y < lO'^ this behavior is reversed. In the third 



Here we show how the initial conditions given in 
Eqs. (fT6|) are determined. We start from the equations 
pi)) for the coefficients Ofc for large k 



2U sin(fc 0^ 



0, 



(Al) 



noting that Uk — C/ for large k (0+ denotes a positive 
infinitesimal increment). This equation is easily solved 
using the initial conditions at /c = 00 given in Eq. (jl5p 
with the result 

akjj = —U ^1 Si(fc 0^)^ Sjj — V6ji + fJ.Sjj, 

ak,j{j±i) = 0. (A2) 

Here, Si(z) is the sine integral function. The impurity 
with strength V is placed at site j = 1. This result can 
now be used in order to solve the equation for a^fl for 
large k, 

L^^-V sin(fc 0+) LU MkO+) f 2 
flfc.o = ; r 1 - -Si(fc 0^) 

TT k TT k \ TT 

(A3) 

where the expansion ln(l + x) ~ x for x < 1 has been 
used. Integration of Eq. (jA3p yields 



y L /U 



where the limit 0+ — > has been performed. 



(A4) 
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APPENDIX B: NUMERICAL SOLUTION OF THE 
FLOW EQUATIONS FOR THE SELF ENERGIES 
AND THE GROUND STATE ENERGY 

The solution of the flow equations for the self-energies 
and ground state energy requires the calculation of 
the inverse and determinant of (potentially) large cyclic 
tridiagonal matrices at each step of the integration of 
a set of ordinary differential equations. Of course, this 
could be achieved with standard library routines. How- 
ever, due to the special form of the Eqs. (fT4|) we only need 
to determine the cyclic tridiagonal part of these matri- 
ces. Therefore, in the following we will develop meth- 
ods adapted to our special needs. As a consequence we 
achieve a significant speed-up of the numerical calcula- 
tion as well as considerable reduction of the computer 
memory requirements. Such improvements enable us to 
treat large system sizes. 



1. Inversion of cyclic tridiagonal matrices 

For matrices M which can be represented as the sum 
of a tridiagonal matrix T and a direct product of two 
vectors u and v. 



M = T + u®v, 



(Bl) 



the inverse can be obtained using the Sherman-Morrison 
formula [53 



z ® w 

1 + V • z 



with 



z^T-^u, w=(T"i)'^v. 
A general cyclic tridiagonal matrix 



(B2) 
(B3) 



M 



/ ai 6i 

■•• ■•■ 







: c„_i a„„i bn^i 
\bn ■■ ■ c„ an J 



(B4) 



can be written in the form (|Bip with the tridiagonal part 
T given by 



f ai bl 

C2 a2 62 

'•• ■•. 



\ 





(B5) 



: c„_i a„_i bn~i 

V • ■ • c„ a„ y 

with ai = 2ai, a„ = a„ -f- ci6„/ai and a; 
1 < i < n. The vectors u and v are defined by 



Ui for 



u = (-fli, 0, • • • , 0, 6„) , V = (1, 0, • • • , 0, -ci/ai) 



(B6) 



For the inversion of the tridiagonal part T, we use 
a method inspired by Andergassen et al. [s^, which is 
based on Ref. [5§]. However, in our case the tridiagonal 
matrix is not complex symmetric, therefore, we have to 
use an LU decomposition of T instead of the LDU de- 
composition employed in Ref. [53|. Details are described 
in the following subsection. 



2. Implementation 



According to Eq. (jB2|l the inversion of a general cyclic 
tridiagonal matrix requires two essential steps: (i) the in- 
version of a tridiagonal matrix and (ii) the determination 
of the vectors z and w. 

(i) A tridiagonal matrix can be represented as a prod- 
uct of two matrices T = LU — UL, where U, U are 
upper diagonal and L, L are lower diagonal matrices. 
From the LU decomposition 



= LU 

C2 d2 



\ / I Ui 

1 U2 



Cn-1 dn-i 

Cn dji / 



(B7) 
\ 



1 Un-l 

1 / 



with 



di=ai, di ^ ai ~ Ui^id, Ui^bi/di (B8) 

one easily finds for 1 < i < j < n the relations 

(T-i)„,„ = l/d„, (T-i),j =-u,(T-i),+i,j. (B9) 

Similarly, from the UL decomposition of T, 

T = UL= (BIO) 
f 1 ui \ / di \ 



1 U2 



1 Un-l 

1 / 



C2 d2 



Cn dji / 



with 



dn = an di ^ a-i - UiCi+i Ui^bi/di+i, (Bll) 
one obtains for I < i < j < n the relations 



(T-i)i,i = 1/di, (T-\, = -{i,_i(T- 



(B12) 



Combining Eqs. (lB9| and ((BT2)) yields a recursion re- 
lation for the diagonal elements of 



(T ^)i+ii+i — — (T ^)ji 

Ui 



(B13) 
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with the initial condition given in Eq. ((BT2I) . From the 
diagonal elements one then generates the upper diagonal 
as well as the upper corner element using the relations 
given above. Alternatively, the UL decomposition yields 
the relation 

(T-i),., =-^(T-i),_i,j forl<j<i<n, (B14) 

which is used to generate the lower diagonal and corner 
element, respectively. 

(a) The vector z is obtained from Eq. (IBSp . Tz = 
LUz = u, as follows: first we solve recursively Lz' = u 
for z' 

z[ — Ui/di, z[ = {ui—Ciz[_^)di for 1 < i < n, (B15) 
and then obtain z recursively from Uz = z' 

Zn = z^, Zi = z,- - MiZi+i for 1 < i < n. (B16) 
In a similar way one determines w from T"^w = v. 







with 60 = ^ri- For large systems, the products in this 
formula can easily underflow or overflow. This must be 
carefully controlled appropriately. 



3. Determinant of a cyclic tridiagonal matrix 



The determinant of a cyclic tridiagonal matrix is cal- 
culated using a formula given in Ref. (60| 



/ai 61 

C2 0.1 &2 



det 











Cn-1 fln-l 

C7) TTTT 



(B17) 
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